function [cb1_L,cb1_R,cb2_L,cb2_R] = func_SUZ(Y,X,T,B,alpha,K_iter,t_grid)

% Parameters
n=length(Y);
p=size(X,2);
tn=length(t_grid);
xi=1+norminv(rand(n,B))/sqrt(2)+(norminv(rand(n,B)).^2-1)/2;
xi=xi./repmat(sum(xi,1),[n,1]);    

% Calculate h_1, h_2, lambda_global, lambda_local
h_1=1.06*std(T)*n^(-1/5); 
h_2_grid=(1+.2*(-1:.1:1))*h_1;
MSE_list = NaN(size(h_2_grid));
for MSE_index = 1:size(MSE_list,2)
    MSE_list(MSE_index) = func_MSE(T,Y,h_2_grid(MSE_index));
end
h_2=n^(-1/10)*1.095*min(h_2_grid(min(MSE_list)==MSE_list));
lambda_global=1.1*norminv(1-(1/log(n))/max(n*h_1,p))*sqrt(n);  
lambda_local=sqrt(log(log(n*h_1)))*sqrt(log(max(n*h_1,p))*n*h_1); 

% Estimate f(x) and v(x)
f_est_T=NaN(n,n);
nu_est_T=NaN(n,n);
ui_T=(repmat(T,[1,n])-repmat(T',[n,1]));
W=max(eps,normpdf(ui_T/h_1)/h_1);
for t=1:n
    temp1=func_Logit(T<=(T(t)+h_1),X,n,p,lambda_global,K_iter);
    temp2=func_Logit(T<=(T(t)-h_1),X,n,p,lambda_global,K_iter);
    f_est_T(:,t)=max(eps,abs(temp1-temp2)/(2*h_1));
    [nu_est_T(:,t),~]=func_LocalReg(Y,X,n,p,h_1,W(:,t),lambda_local,K_iter);
end

% Estimate mu(Ti) and mu_b(Ti)
W=max(eps,normpdf(ui_T/h_2)/h_2);
mu_est_T=NaN(n,1);
mu_boot=NaN(n,B);
for t=1:n    
    Pi_n=(Y-nu_est_T(:,t)).*(W(:,t)./f_est_T(:,t))+nu_est_T(:,t);    
    mu_est_T(t)=mean(Pi_n,1);     
    mu_boot(t,:)=Pi_n'*xi;
end

% Estimate the standard error 
ui_t=(repmat(T,[1,tn])-repmat(t_grid,[n,1]));
f_est_t=NaN(n,tn);
nu_est_t=NaN(n,tn);
W=max(eps,normpdf(ui_t/h_1)/h_1);
for t1=1:tn
    temp1=func_Logit(T<=(t_grid(t1)+h_1),X,n,p,lambda_global,K_iter);
    temp2=func_Logit(T<=(t_grid(t1)-h_1),X,n,p,lambda_global,K_iter);
    f_est_t(:,t1)=max(eps,abs(temp1-temp2)/(2*h_1));
    [nu_est_t(:,t1),~]=func_LocalReg(Y,X,n,p,h_1,W(:,t1),lambda_local,K_iter);
end
sig_est=NaN(tn,1);
Kbar=(ui_t.*exp(-(ui_t.^2)/4))./(4*sqrt(pi));
for t1=1:tn
    Pi_n=((Y-nu_est_t(:,t1))./f_est_t(:,t1).*(Kbar(:,t1))).^2;    
    sig_est(t1)=sqrt(mean(Pi_n,1)/h_2);     
end

% Calculate beta1(t) and beta1_b(t) 
W=max(eps,normpdf(ui_t/h_2)/h_2);
b1_est=NaN(tn,1);
b1_boot=NaN(tn,B);
b1_diff=NaN(tn,B);
b2_diff=NaN(tn,B);
for t1=1:tn
    beta_t=lscov([ones(n,1),T-t_grid(t1)],mu_est_T,W(:,t1));
    b1_est(t1)=beta_t(2);
    for b=1:B
        beta_b=lscov([xi(:,b),xi(:,b).*(T-t_grid(t1))],xi(:,b).*mu_boot(:,b),W(:,t1));
        b1_boot(t1,b)=beta_b(2);
        b1_diff(t1,b)=sqrt(n*h_2^3)*(b1_boot(t1,b)-b1_est(t1))/sig_est(t1);
        b2_diff(t1,b)=sqrt(n*h_2^3)*abs(b1_boot(t1,b)-b1_est(t1))/sig_est(t1);
    end      
end

% Calcualate the confidence band
b1_max=max(b1_diff);
cv1=quantile(b1_max,1-alpha);
cb1_L=b1_est-cv1*sig_est/sqrt(n*h_2^3);
cb1_R=b1_est+cv1*sig_est/sqrt(n*h_2^3);

b2_max=max(b2_diff);
cv2=quantile(b2_max,1-alpha);
cb2_L=b1_est-cv2*sig_est/sqrt(n*h_2^3);
cb2_R=b1_est+cv2*sig_est/sqrt(n*h_2^3);

end

